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Abstract 

To study the ballistic transport of charge carriers in nano-structured quantum devices, a highly 
efficient numerical technique is developed, which provides continuous transmission spectra for ar- 
bitrarily complex potential geometries in two dimensions. We apply the proposed method to 
single and double barrier structures and compare the results with those obtained using standard 
techniques for computing transmission coefficients. Excellent numerical agreement as well as con- 
siderable computational saving is demonstrated. 
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I. INTRODUCTION 



A wide range of potential applications for nano-scale electronics in the quantum regime, 
together with a rapid increase in computing power, have generated much interest in the nu- 
merical analysis of nano-structured devices as a viable means for studying their properties 



. Examples of such applications based on the transport properties of ballistic 
charge- carriers include quantum wires 
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nelling diodes and band filters LUA LL^U quar 
„ .Witches Q.,uan.u..en.o.lyQ 



quantum resistors |10], resonant tun- 



uantum transistors and stub tuners 
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gates 
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as well as qubits and quantum logic 



2l|. 



Much insight into the behavior of such devices can be gained by studying the energy 
dependance of transmission through various nano-structures for propagating charge carri- 
ers. The transmission coefficient of a propagating wave function can be determine d by 
the application of time-independent (static) or time-dependent (dynamic) techniques |22|. 
The time-independent approach involves solving the time-independent Schrodinger equa- 
tion, commonly carried out using one of three methods: (1) Mode matching method which 
is based on splitting the nano-structures into different regions with known analytical solu- 
tion. The total wave function is then derived by matching the solutions at the boundary 
between re gion s. This method has been widely used to study quantum devices with simple 



geometries 
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y, 25, Q- 

In principle, any arbitrary structure can be divided into many 
segments, each of which is approximated as a rectangular potential well with a finite width. 
For complex geometries however, this can result in a very large matrix to be inverted and the 
computation can become unstable. (2) Finite element method 23| where using an irregular 
numerical mesh, the potential and the wave function are discretized by choosing a suitable 
set of basis functions to approximate each grid point. However, due to the approximate 
nature of the elemental solutions, large grids are still required for accurately representing 
the wave function particularly when dealing with excited or continuous energy states. This 
often results in an excessive computational cost as well as the introduction of high frequency 
noise, which can considerably distort the phase of the wave function. (3) Green's function 



method 



m 



29[ where the Schrodinger's equation is expressed and solved using the common 



Green's function techniques. Although in principle capable of providing solutions to arbi- 
trarily complex problems independently of the potential boundary structures, the Green's 
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function method is most efficient when used in conjunction with simphfying assumptions 
such as the decouphng of the x and y variables in the two-dimensional Schrodinger equation 
or by not requiring information about a global wave function. In particular, the later case 
leads to an efficient algorithm for calculating the transmission coefficients known as the re- 
cursive Green's function method j^, where it is necessary to compute the Green's function 
only for some sections of the structure. The calculation of the wave function for a general 
two-dimensional case however relies on computing the generalized Green's function for the 
entire region, which is computationally expensive. 

As well as their limited utility in constructing global wave functions for complex nano- 
structure geometries efficiently, the time-independent methods are also unable to provide 
information on the transient behavior of the system under study. This can be addressed 
by solving the time-dependent Schrodinger equation where transmission coefficients are ob- 
tained by numerically propagating an initial wave function through a given potential geom- 
etry and then summing over the transmitted parts. Methods for solving the time-dependent 
Schrodinger equation include: (1) The finite difference method Isil. Is^. which uses an nth 
order finite difference approximation to expand the time-evolution operator. This method 
however suffers from a low convergence rate as well as excessive truncation errors over mod- 
estiv ,o.g propagation t„„e. (2) Sp„t ope.ato. .nethod 0. whe.e the ti.e-evolution 
operator is approximated as the product of three diagonal operators which can be readily 
solved. This scheme however neglects the commutator between these operators and thus 
introduces error in both e nerg v and phase of the wave function. (3) The time-dependent 
Green's functions method |34], which employs the standard Green's function techniques. 
Due to its reliance on summations over a large number of energy-eigenstates however, this 
method too becomes computationally prohibitive when dealing with complex potentials. (4) 
Chebyshev-Fourier propagation scherne, which is based on expressing the evolution operator 
using the Chebyshev approximation This method proves to be numerically robust as 
well as efficient for propagations over arbitrary potentials. 

Nonetheless, when approached naively, calculating the transmission coefficients for a wide 
range of energies using the Chebyshev-Fourier scheme could still prove computationally 
expensive, amounting to repeated propagations of an initial wave function with various mean 
energies in the given range. A novel method for computing the transmission coefficients was 
developed by Yin and Wang [2^, which arrives at the entire transmission spectrum for a 
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given range of energies after only one propagation of an initial wave function with a broad 
energy distribution. Although this method is shown to be highly efficient while maintaining 
the capacity to deal with arbitrary potentials in one dimension, its use in two dimensions 



is limited to geometries wit 
decoupled x and y variables 



1 infinite boundary conditions in the second dimension and 
36, 



making it possible to treat 



In this paper we extend the work of Yin and Wang 
arbitrarily complex potentials in two dimensions, without imposing any limitations on the 
dimensional parameters of the system, such as decoupling in the x and y directions. We apply 
this scheme to single and double hyperbolic barrier structures and demonstrate excellent 
numerical accuracy while maintaining a high degree of computational efficiency, compared 
to standard techniques for obtaining individual transmission coefficients. 

II. THEORY 

We begin by considering the two-dimensional Schrodinger equation 

d 

iH—i>{x, y, t) = Hipix, y, t), (1) 

where the Hamiltonian H is given by 

H = -^V' + Vix,y,t) (2) 
2m* 

with m* being the effective mass of the propagating charge carrier within the semiconducting 
lattice. If the potential function is time-independent i.e V{x,y,t) = V{x,y), then a formal 
solution to this equation is 

t ^ 

ij{x, y, t) = exp{--Ht)ij{x, y, 0) (3) 

where t is the propagation time, and exp{—j^Ht) is the time evolution propagator, commonly 
denoted by U. The Chebyshev-Fourier scheme as detailed in approximates f/ by a 

Chebyshev polynomial expansion 

ipix, y, t) = exp[-i{£rnax + £min)t] a„(a)0„(-7^)V^(x, y, 0) (4) 

n=0 

where Smin and Smax are the upper and lower bounds on the energies sampled by the wave 
packet, an{a) = 2Jn{a) except for 00(0;) = Jo{a), Jn{a) are the Bessel functions of the first 



kind, (pn are the Chebyshev polynomials, and TV is the number of terms in the Chebyshev 
expansion. To ensure convergence, the Hamiltonian needs to be normalized as 

7; 7^ £max ^min] • (5) 

The action of the Laplacian operator V on the wave functions is carried out using a Fourier 
transformation technique It is important to note that the rapid convergence and the 

ligh n umerical accuracy of this propagation scheme are well established in the literature 
4, 35, 38, 39, 4^, where it has been used in the study of various quantum nano-structures. 

In the simulations we start with a localized wave packet as the initial wave function 
ipi = ip{x, y,t = 0) and then use the Chebyshev-Fourier scheme as a means to propagate this 
wave packet in time t. The final wave function ipf = ilj{x,y,t = t) is partly transmitted and 
partly reflected due to the potential barrier V{x,y) (Fig. Q). To evaluate the transmission 
coefficient T{S) for any given potential barrier, one can setup ipi with an incident mean 
energy S and a very small energy uncertainty, propagate ipi in time until it is completely 
scattered by the barrier, and then collect the transmitted parts oi ipf. In this way, obtain- 
ing a transmission spectrum entails multiple propagations for different incident energies as 
depicted in Fig. |21 Given the computational cost associated with any time-dependant prop- 
agation scheme however, the use of this direct method for obtaining a continues transmission 
spectrum is clearly inefficient. 

In Yiu and Wang's approach 3^, a wave packet with a small spatial uncertainty (cor- 
responding to a large momentum spread) is used to obtain the entire transmission carve 
after only one propagation. In other words, one can compute the transmission coefficients 
by simply transforming the transmitted wave packet to momentum space and dividing by 
the original momentum space wave function. Noting that in one dimension, the transmitted 
and reflected wave packets can be expressed as momentum space components of ipf with 
positive momentum {p > 0) and negative momentum {p < 0) respectively, one can then 
write 

fi(rtJ^V(P<°'l\ (6) 

where T{p) and R{p) are the transmission and reflection coefficients at momentum p. This 
scheme however, is strictly limited to scattering in one dimension. In order to show this we 
present a more generalized approach to the problem, which enables us to extend the current 
scheme to two dimensions. 
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First we note that in principal ip{x^y,t) may be expressed as a linear superposition of 
energy eigenfunctions (f)E{x,y). We can therefore rewrite Eq. Q as 



^lJ{x,y,t) = e-^^*/"^Q(i?)0s(x,2/) (7) 

= ^e~'^'/''Q{E)M^,y) (8) 

= ^CfiE)M^,y) , (9) 

where Ci{E) and Cf{E) denote the initial and final complex coefficients associated with 

each energy eigenfunction. We note that Ci{E) and Cf{E) only differ by a phase factor. 
Expressing Cf{E) as {E\ipf) we can expand Cf{E) using the position basis set 

Cf{E) = Jj {E\x,y){x,y\iJi) dxdy (10) 

\cf{E)f = [[ {E\x,y){x,y\^f){i;f\E) dxdy (11) 



In the position space however the final wave function assumes two components: the trans- 
mitted part lying in an area T, and the refiected part lying in an area R. Noting that these 
two do not overlap we can write 



\cjiE)f = jj^{E\x,y){xM^f)hPf\E)dxdy (12) 
+ fj {E\x,y){x,y\iPf){^f\E) dxdy (13) 



J JR 

= \cT{E)f + \cR{E)f (14) 

where ct{E) and cr{E) are the transmission and refiection coefficients of each energy eigen- 
function. We therefore conclude that 

UE)\' = \cf{E)\' = \cT{Et + \cr{E)\\ (15) 

In other words the amplitude of each energy eigenfunction will independently conserve 
throughout the propagation and must therefore have its own transmitted and reflected parts. 

In order to evaluate |q(-E')|^, |ct'(-E')|^ and |c/j(i?)|^, we consider the case where the initial 
and the scattered wave functions are both sufficiently distant from the potential barrier, 
such that the interaction between the wave function and the potential barrier is negligible, 

6 



i.e. V ^0. As a result we can ignore the contribution of the potential term V{x, y) 
in the Hamiltonian since H p^/2m*, meaning that energy is directly proportional to the 
square of the momentum operator. 

Now representing the initial wave functions in the momentum space, denoted by 
i'iiPxyPy) = 1 ^^Pii^ (Px^ + Pyy))4'i{x , y)dxdy , and noting that the total momentum 

p"^ = pI + pf/y "we have 

\(^^iE)f = ^'^\^iiP-^Py)f^ ioTpl+pl = 2m*E. (16) 

In this way computing \ci{E)f amounts to simply summing over those momentum com- 
ponents which lie on a circle with radius V2m*E as depicted in Fig. El Furthermore, by 
defining the motion of the initial wave packet to be in the positive x direction (see Fig. 

transmission and reflection can also be conveniently characterized using the momentum 
space representation of the final wave function ipfiPx^Py)- Referring to the momentum circle 
described above, \cT{E)f and \cR{E)f correspond to semicircles with Px > and px < 



respectively, i.e. 



\cT{E)\'=^^\MPx,Py)\\ hTpl + pl = 2m*E, (17) 

Px>0 Py 

\cRiE)\'=^^\^jiPx,Py)\\ ior pl+pl = 2m*E. (18) 

Px<0 py 

In this way, the transmission coefficient T{E) can be written as 

T(E) = . (19) 

Therefore, in principle, all the transmission coefficients can be calculated after propagating 
only a single initial wave packet ipi for any choice of In this work we have used an initial 
wave packet which has a small momentum spread in the y direction, but a broad momentum 
distribution in the x direction, as depicted in Fig. |3J The scattering of this customized wave 
packet, which we refer to as the broad- energy wave packet, provides us with an efficient means 
for computing the continuous transmission spectrum over a wide range of incident energies 
for an arbitrary potential barrier in two dimensions. We term this scheme the momentum 
space method. 
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III. RESULTS 



In this section we demonstrate the use of the momentum space technique to obtain a series 
of continuous transmission spectra for the scattering of balhstic electrons by the single and 
double barrier potentials (Fig.E)), positioned at various angles a with respect to the incident 
electron (Fig. P). We then compare the results with transmission coefficients obtained using 
the direct method (using multiple wave packets with small momentum uncertainty). 

The barriers are constructed using hyperbolic functions in the x direction: V{x) = 
Vo/ cosh^(^) for the single barrier and V{x) = Vo/ cosh^(^) + Vq/ cosh^(^) for the double 
barrier where the barrier height Vq = 13.6 eV, width parameter a = 0.21 nm (giving a half 
width of 0.19 nm), and the barrier separation 2d = 1.58 nm. Simulations are carried out 
using three orientation angles a = 90°, 62°, and 52°. 

The incident ballistic electrons are modelled using a Gaussian wave packet 

V'(x, y, 0) = expii-ix - Lf /2al - /2al) + i{xp, + ypy)/h) , (20) 

where ax = h/\/2^px and ay = h/ \/2Apy are the wave packet's standard deviation in the 
X and y dimensions, ^px and Ap^ are momentum uncertainty components, p^ 
and = are the mean momentum components of the wave packet with mean energy 
£. In the simulations we consider electron wave packets with effective mass m* = 0.06me 
and £ in the range of [0.1,3.5]Vo. The narrow-energy wave packets (used in the direct 
method) have a spacial distribution given by ax = ay = 1 .'b nm corresponding to /S.E in 
the range of [0.02, 0.1]Vo, and L = 45.0 nm. The broad-energy wave packets (used in the 
momentum space method) are characterized by ax = 0.75 nm and ay = 7.5 nm, mean energy 
S = IVq, and L = 26.5 nm. Wave packets are propagated for a time t = rriD/px where 
D ^ 3L is the propagation distance. The choice of parameters L and D ensures that both 
the initial and final wave functions are sufficiently away from the potential barriers (i.e. 



V 



ipj/i^ H ipj < 10""^)). Fig. ini and m illustrate the position space evolution of the 
broad-energy wave packet for the double barrier potentials with a = 90° and 52°. 

In order to examine the accuracy of the proposed momentum space method, we show 
firstly that the transmission curves are independent of the choice of S for the broad-energy 
wave packets. Fig. |H1 presents a transmission curve for the double barrier potential with 
a = 52°, constructed by overlapping multiple segments from transmission curves obtained 
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using initial wave packet with £/Vq = 0.2,0.6,1.0,1.4,1.8,2.2,2.6, and 3.0 respectively. 
An excellent match between these segments producing a uniform transmission spectrum is 
demonstrated. 

Secondly, we compare the continuous transmission curves obtained using the momen- 
tum space method, against the discrete transmission coefficients obtained using the direct 
method. Fig. El and [TUl depict the results for the single and double barrier potentials re- 
spectively, with three orientation angles a = 52°, 62° and 90°. A close match between each 
transmission curve and its corresponding transmission coefficients (points on each curve) 
demonstrates an excellent agreement between the two methods. 

In order to examine the computational efficiency of the momentum space method, we 
ffist note that the accuracy and convergence of the Chebyshev- Fourier propagation scheme 
primarily depends on the choice of Smax in Eq- which determines the dimensions of the 
numerical grid due to Nyquist's sampling theorem. This accuracy improves asymptotically 
with higher values of Smax while making computation more costly. In the simulations Smax 
is estimated based on the maximum potential energy of the system Vq, the mean momentum 
of the initial wave packet, and its maximum momentum range. We then define an accuracy 
factor f] to be the ratio between the actual value of Smax used in the simulation and the ini- 
tial estimation. All results from the momentum space method were obtained using f3 = 1.0, 
which was found to provide good accuracy while maintaining a reasonable computational 
speed. The direct method however, at times, required a higher (3 value for convergence. 
In the case of the double barrier with a = 90° for instance, increasing /3 from 1.0 to 1.2, 
produces no discernable improvements in the resulting transmission curve using the momen- 
tum space method (Fig. fTT| . Using the direct method however, setting /5 = 1.0 leads to 
significant discrepancies between the transmission points and the transmission curve (Fig. 
IT^ . Conversely, setting (3 = 2.0 results in an excellent agreement (Fig. fT^ . but severely 
increases the cpu time requirement. 

The efficiency of the momentum space scheme becomes particularly apparent when con- 
sidering the existence of narrow transmission peaks in the spectrum (such as the one at 
S = OAVq for the double barrier potential with a = 52°). In the above simulations, obtain- 
ing a continuous transmission curve using the momentum space method was 5 times faster 
than the application of the direct method, when requiring only 20 points to represent the 
same curve. In the absence of prior information about the topography of the transmission 
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spectrum however, using the direct method necessitates probing the spectrum with a much 
higher resolution, in order to detect possible narrow transmission peaks. As such, apply- 
ing the direct method can become prohibitively inefficient, rendering the use of momentum 

space scheme vastly advantageous. 

IV. CONCLUSION 

We presented an efficient scheme for computing the continuous transmission spectrum 
of charge carriers scattered by an arbitrary potential geometry in two dimensions. We 
applied this scheme to single and double barrier potentials at various angles of incidence, and 
demonstrated excellent precision as well as substantial computational saving as compared 
to a standard method commonly used for calculating transmission coefficients. 
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V. FIGURES 
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FIG. 1: Position space: the incident particle approaches the scattering potential barrier from the 
left in the x direction. The angle a describes the orientation of the barrier with respect to the 
incident particle. 
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FIG. 2: Momentum space: transmission and reflection are defined as wave components with px > 
and Px < respectively. Components il^{E) correspond to rings of radius r = V2m*E. 
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FIG. 3: In the momentum space incident charge carriers are modelled using initial wave packets 
with a narrow momentum spread, positioned at m*£ and Py = 0, where £ is the mean 

energy of the wave packet. A complete transmission curve can be constructed by computing the 
transmission coefficient of every wave packet with £ in the given range. 




FIG. 4: A broad-energy wave packet has a large momentum spread in the x direction, corresponding 
to the energy range of the transmission spectrum. In the y direction however, it has the same spread 
as the single charge carrier wave packets. 
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FIG. 7: Snapshots of the broad-energy wave packet scattering from the double barrier potential 
with Q = 52°. 
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FIG. 8: Continuous transmission spectrum for the double barrier potential with a = 52°. Eight 
broad-energy wave packets with mean energies £ /Vq = 0.2,0.6,1.0,1.4,1.8,2.2,2.6, and 3.0 have 
been used to obtain the corresponding segments of the final curve. The agreement between these 
overlapping segments demonstrates that the transmission curve is independent of the choice of £. 




E/VO 

FIG. 9: Continuous transmission spectra for the single barrier potential with orientation angles 
a = 52° (solid), 62° (dashed) and 90° (dotted). Points on each curve represent the corresponding 
transmission coefficients calculated using the direct method. 
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FIG. 10: Continuous transmission spectra for the double barrier potential with orientation angles 
a = 52° (solid), 62° (dashed) and 90° (dotted). Points on each curve represent the corresponding 
transmission coefficients calculated using the direct method. 
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FIG. 11: Transmission curves for the double barrier potential with a = 52°, obtained using the 
accuracy factor /3 = 1.0 (solid) and /3 = 1.2 (dotted). The agreement between the two curves 
demonstrates that good convergence has already been achieved at /? = 1.0. 
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FIG. 12: Transmission curve for the double barrier potential with a = 52° vs. direct method 

transmission points. Both results have been obtained using the accuracy factor (3 = 1.0. Differences 
between the transmission curve and its corresponding transmission points suggest that the direct 
method does not achieve complete convergence at /? = 1.0. 
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FIG. 13: Transmission curve for the double barrier potential with a = 52° vs. direct method 
transmission points. By increasing the accuracy factor /? of the direct method from 1.0 to 2.0, we 

arrive at an agreement between the transmission curve and its corresponding transmission points. 
This suggests that the direct method requires /3 > 1.0 for complete convergence. 
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